library(dplyr)
library(sf)
library(readr)
library(mapview)
library(Rspatialworkshop)
library(ggplot2)
library(leaflet)
library(httr)
Use read_csv fromreadr package to load a
.csv file I keep as part of a spatial data package that you can install
with:
library(devtools)
install_github("mhweber/Rspatialworkshop")
library(Rspatialworkshop)
Then we:
gages <- read_csv(system.file("extdata/Gages_flowdata.csv", package = "Rspatialworkshop"),show_col_types = FALSE) %>% dplyr::select(ID=SOURCE_FEA,LON_SITE,LAT_SITE, MeanFlowCFS=AVE, DrainAreaSqMiles=DA_SQ_MILE) %>%
dplyr::filter(DrainAreaSqMiles > 500)
glimpse(gages)
Rows: 539
Columns: 5
$ ID <dbl> 14361500, 14378000, 14335000, 1…
$ LON_SITE <dbl> -123.3178, -123.8123, -122.5956…
$ LAT_SITE <dbl> 42.43040, 42.37900, 42.69985, 4…
$ MeanFlowCFS <dbl> 3402.544, 2317.657, 1796.864, 7…
$ DrainAreaSqMiles <dbl> 2459, 665, 650, 698, 1610, 670,…
st_as_sf from the sf package to
promote the data to a simple feature collection by:
epsg codeepsg code is for unprojected NAD83gages <- st_as_sf(gages, coords = c("LON_SITE", "LAT_SITE"), crs = 4269, remove = FALSE)
ggplot() + geom_sf(data=gages)
mapviewmapview(gages)
m <- mapview(gages, map.types = c("Esri.WorldShadedRelief", "OpenStreetMap"), color = "grey40",cex = "MeanFlowCFS")
m
Web Map
Services provide pre-rendered map tiles at different scales and are
useful as background map layers. Here is my original Stack
Overflow question from when I first played with adding WMS layers
with mapview and leaflet
Here we display a background of National Hydrography Dataset stream lines to go with our gages without having to load local files.
This is a super handy way to share your data in an interactive map mashed up with relevant background layers in one file in R, without having to share spatial files.
m <- mapview(gages, cex = "MeanFlowCFS")
# Define the NHD WMS service
wms_nhd <- "https://hydro.nationalmap.gov:443/arcgis/services/nhd/MapServer/WmsServer?"
m@map <- m@map %>%
addWMSTiles(group = 'USGS HydroCache',
wms_nhd,layers = 6,
options = WMSTileOptions(format = "image/png", transparent = TRUE),
attribution = "") %>% mapview:::mapViewLayersControl(names = "NHD")
m
We can see other services available on the National Map here - let’s try transportation
# Define the WMS service
m <- mapview(gages, cex = "MeanFlowCFS")
wms_trails <- "https://carto.nationalmap.gov:443/arcgis/services/transportation/MapServer/WmsServer?"
m@map <- m@map %>%
addWMSTiles(group = 'Trails',
wms_trails,layers = 3,
options = WMSTileOptions(format = "image/png", transparent = TRUE),
attribution = "") %>%
addWMSTiles(group = 'Trail Labels',
wms_trails,layers = 15,
options = WMSTileOptions(format = "image/png", transparent = TRUE),
attribution = "") %>% mapview:::mapViewLayersControl(names = c("Trails","Trail Names"))
m
Just to show scope of different things we can do linking to REST Services as well - here are links to a couple examples- spatial wfs services and accessing arcgis rest services
And we can also use REST services for the same WMS layers we displayed in maps above! See National Map Services page
We can get a listing of all ESRI REST services here
url <- parse_url("https://services.arcgis.com/P3ePLMYs2RVChkJx/arcgis/rest/services")
url$path <- paste(url$path, "USA_Water_Bodies/FeatureServer/0/query", sep = "/")
url$query <- list(where = "STATE = 'OR'",
outFields = "*",
returnGeometry = "true",
f = "geojson")
request <- build_url(url)
wb <- read_sf(request)
mapview(gages) + mapview(wb, col.regions='light blue')
This time use the REST service and a query for ‘Pacific Crest National Scenic Trail’ - bonus add a default topography brackground with it
url <- parse_url("https://carto.nationalmap.gov/arcgis/rest/services")
url$path <- paste(url$path, "transportation/MapServer/37/query", sep = "/")
url$query <- list(where = "Name = 'Pacific Crest National Scenic Trail'",
outFields = "*",
returnGeometry = "true",
f = "json")
request <- build_url(url)
PCT <- read_sf(request)
mapview(PCT, color='dark green',map.types = "OpenTopoMap")